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Abstract 

The  task  of  using  molecular  dynamics  to  determine  the  shear  viscosity  and  the  thermal 
conductivity  of  a model  fluid  is  examined  in  detail  for  the  Lennard-Jones  fluid.  The 
equilibrium  time  correlation  function  approach  is  compared  with  three  nonequilibrium 
methods,  Poiseuille  Flow,  Momentum  Impulse  Relaxation,  and  Reversed  Perturbation 
Nonequilibrium  Molecular  Dynamics.  The  conventional  wisdom  is  that  the  equilibrium 
approach  requires  very  long  simulation  times  in  order  to  obtain  statistically  significant 
results.  This  study  finds  that  only  the  Reversed  Perturbation  Nonequilibrium  Molecular 
Dynamics  method  has  the  possibility  of  significantly  reducing  the  simulation  time  required 
to  obtain  “good”  estimates  for  the  thermal  conductivity.  All  of  the  reliable  methods 
examined  can  produce  accurate  estimates  for  the  shear  viscosity  with  shorter  simulation 
times  than  are  needed  for  the  determination  of  the  thermal  conductivity. 

1.  Introduction 

In  this  note,  we  consider  some  schemes  for  estimating  transport  coefficients  of  fluids  that 
do  not  make  use  of  the  Green-Kubo  formulation  of  transport  coefficients  [1].  The  Green- 
Kubo  methods  are  derived  from  statistical  mechanics  and  yield  correct  values  when  prop- 
erly evaluated.  However,  they  require  very  long  run  times  to  obtain  reliable  estimates. 
In  the  alternative  methods  discussed  below,  the  system  is  arranged  so  that  the  proper- 
ties of  the  system  can  be  mapped  onto  a calculable  hydrodynamic  flow  and  the  transport 
coefficient (s)  are  then  obtained  as  parameters  that  make  the  computations  match  the  hy- 
drodynamic solution,  just  as  is  done  experimentally.  The  objective  of  this  study  is  to 
identify  computationally  efficient  methods  for  estimating  transport  coefficients  for  fluids 
that  are  defined  in  terms  of  molecular  structure  and  interactions.  The  methods  examined 
here  are  labeled  Poiseuille  Flow,  Momentum  Impulse  Relaxation,  and  Reversed  Perturba- 
tion Nonequilibrium  Molecular  Dynamics.  In  all  of  the  simulations  discussed  in  this  note, 
the  equations  of  motion  were  integrated  using  the  Beeman  algorithm  [2].  A neighbor  table 
was  used  to  improve  the  efficiency  of  the  simulations  [3]  with  the  result  that  the  simulation 
time  approximately  scales  as  the  number  of  particles. 
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In  order  to  compare  results  for  the  methods  to  be  discussed,  we  list  in  Table  1 some 
Green-Kubo  based  estimates  of  the  transport  coefficients  thermal  conductivity  A and  shear 
viscosity  77  of  the  Lennard- Jones  fluid  for  several  dense  fluid  state  points  [4].  As  usual, 
all  properties  are  listed  in  terms  of  the  Lennard-Jones  potential  parameters  for  energy  (e) 
and  length  (a)  and  in  terms  of  the  mass  m of  the  Lennard-Jones  particle  and  are  therefore 
dimensionless.  The  Lennard-Jones  potential,  </>(r),  has  the  form 

4>(r ) — 4e[(cr/r)12  — (tr/r)6] . (1.1) 

The  time  unit  r for  simulations  is 


r = 


yjmcr2 /e. 


(1.2) 


The  tests  of  the  methods  reported  here  are  all  for  the  Lennard-Jones  fluid  and  the  values 
in  Table  1 are  used  to  evaluate  the  effectiveness  of  the  alternative  methods.  Other  sets 
of  transport  coefficients  for  the  Lennard-Jones  fluid  determined  using  the  Green-Kubo 
method  are  found  in  the  literature  [5,  6,  7,  8,  9].  This  is  not  an  all-inclusive  list  of 
sources.  While  the  Lennard-Jones  fluid  is  an  idealization,  it  has  been  shown  that  when 
the  potential  parameters  are  appropriately  chosen,  the  model  predicts  values  for  the  shear 
viscosity  and  thermal  conductivity  for  the  rare  gases  argon,  krypton,  and  xenon  that  are 
in  good  agreement  with  experimental  values  [10]. 

Table  1.  A set  of  Green-Kubo  based  thermal  conductivity  A and  shear  viscosity  rj  transport 
coefficients  for  the  Lennard-Jones  fluid  at  specified  number  densities  n and  temperatures  T are 
listed  here.  The  estimated  uncertainty  in  the  transport  coefficients  is  reported  to  be  on  the  order 

of  5%  [4]. 


n 

T 

A 

V 

0.85 

0.745 

7.66 

3.28 

0.85 

0.908 

7.18 

2.64 

0.85 

1.033 

7.34 

2.75 

0.85 

1.318 

8.03 

2.56 

0.85 

1.503 

8.19 

2.44 

0.85 

1.976 

8.08 

2.28 

0.85 

2.905 

8.46 

2.31 

0.75 

1.008 

4.72 

1.66 

0.75 

1.330 

5.64 

1.68 

0.75 

2.017 

6.06 

1.51 

0.75 

2.913 

7.02 

1.53 

0.65 

1.588 

4.14 

0.88 

0.65 

1.851 

4.46 

0.99 

0.65 

3.564 

5.30 

1.00 
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The  three  nonequilibrium  methods  are  discussed  in  the  next  three  sections.  In  each  case, 
the  method  is  first  reviewed  based  on  the  existing  literature.  Then  the  way  we  have  imple- 
mented the  method  is  described  and  some  results  for  the  Lennard-Jones  fluid  are  presented. 
The  final  section  contains  a discussion  of  the  merits  of  the  methods  and  some  recommen- 
dations on  which  methods  are  potentially  useful  for  the  simulation  of  more  complex  fluids. 

In  order  to  make  explicit  comparisons  of  the  computational  effort  involved  in  the  different 
methods,  we  first  introduce  the  equilibrium  time  correlation  function  approach  known  as 
the  Green-Kubo  method  and  present  results  for  the  first  state  point  listed  in  Table  1. 

The  Green-Kubo  method  expresses  transport  coefficients  as  time  integrals  of  equilibrium 
time  correlation  functions  of  fluctuating  currents  [11,  12].  For  the  shear  viscosity,  the 
current  is  the  off-diagonal  element  of  the  stress  tensor  Vap, 


N 


^Pa/3  ~ 


ot  ft 
m.ivrVi 


i= 1 


- E'li'T 

j>i 


(1.3) 


The  shear  viscosity  is 

1 f°° 

V = yfcgj1  J < 'Patpity'Pai s{t)  > dt  (1-4) 

where  < . . . > indicates  an  ensemble  average  of  the  enclosed  quantity.  The  corresponding 
current  for  the  thermal  conductivity,  the  heat  current  Ja , takes  the  form  in  the  micro- 
canonical  ensemble  [12], 


Ja  = ^2  \vf  [\rnvf  + 4>{rij)]  + 1 • Vi)rf . 

i=  1 j^i  j^i 

The  thermal  conductivity  is 


(1.5) 


< Ja(t)Ja(t)  > dt.  (1.6) 

It  is  important  to  realize  that  the  heat  current  expression  is  ensemble  dependent  [12]. 

The  shear  viscosity  and  thermal  conductivity  have  been  determined  for  a 500  Lennard- 
Jones  particle  system  in  the  microcanonical  ensemble.  That  is  for  constant  number  of 
particles  N,  constant  system  volume  V,  and  constant  energy  E.  The  total  momentum  of 
the  system  is  set  to  zero.  A total  of  106  time  steps  of  duration  0.01  r were  needed  to 
obtain  satisfactory  convergence  of  the  current-current  time  correlations  to  zero  for  long 
times.  This  amounts  to  an  overall  time  interval  of  104r.  A total  of  3985  time  origins 
separated  by  2.5r  were  generated  so  that  the  time  correlation  functions  are  averaged  over 
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3985  independent  samples.  The  influence  of  system  size  on  the  transport  coefficients  is  less 
than  the  uncertainty  in  the  estimates  for  N>500  [13,  14,  15,  16].  In  addition  to  the  explicit 
results  for  the  state  with  T = 0.76  and  n = 0.85  that  are  presented  below,  the  other  states 
listed  in  Table  1 have  been  run  in  order  to  estimate  the  uncertainties  in  the  transport 
coefficients.  A rough  estimate  of  the  uncertainty  in  these  coefficients  can  be  obtained 
by  determining  the  variance  of  the  values  of  the  integrals  of  the  three  time  correlation 
functions  [17].  In  general  the  uncertainties  in  the  shear  viscosity  are  on  the  order  of  5% 
and  those  in  the  thermal  conductivity  are  on  the  order  of  10%. 

Shear  viscosity 

The  time  correlation  functions  were  generated  for  the  xy,  xz,  and  yz  components  of  the 
stress  tensor  and  are  shown  in  Fig.  1.1  as  Vaf3(t)  scaled  so  that 


V = 


Va0(t)dt. 


(1.7) 


P“P(t) 

ap:  xy,  xz,  yz 


t/x 

Fig.  1.1.  The  normalized  time  correlation  functions  that  appear  in  eq.  1.5  are  shown  as  functions 
of  time. 

Fig.  1.2  shows  the  average  of  the  three  time  correlation  functions  and  the  integral  of  the 
averaged  correlation  function.  The  “long  time”  value  of  the  integral  is  3.38,  in  reasonable 
agreement  with  the  value  listed  in  Table  1.  The  standard  uncertainty  for  this  state  is 
±0.2. 
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Fig.  1.2.  The  averaged  time  correlation  functions  (solid  line)  and  the  time  integral  for  the  shear 
viscosity  (dotted  line)  are  shown  here. 

Thermal  conductivity 

The  current-current  correlation  functions,  Ja(t),  scaled  so  that 

/•OO 

A = / Ja(t)dt  (1.8) 

Jo 

are  shown  in  Fig.  1.3  for  the  x,  y,  and  z components  of  the  heat  current. 


J„(t) 


t/x 


Fig.  1.3.  The  averaged  time  correlation  functions  for  the  heat  current  are  shown  here. 

Fig.  1.4  shows  the  average  of  the  three  time  correlation  functions  and  the  integral  of  the 
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Fig.  1.4.  The  averaged  time  correlation  functions  (solid  line)  and  the  time  integral  for  the  thermal 
conductivity  (dotted  line)  are  shown  here. 


averaged  correlation  function.  The  “long  time”  value  of  the  integral  is  7.50,  in  reasonable 
agreement  with  the  value  listed  in  Table  1.  The  standard  uncertainty  in  this  value  is  ±0.8. 

Note  that  the  long  time  limit  for  the  time  correlation  functions  for  this  state  point  is 
reached  when  t < 1.5r  so  that  the  individual  samples  of  the  correlation  functions  are 
independent. 

2.  Poiseuille  Flow 

Poiseuille  flow  refers  to  the  flow  of  a fluid,  subject  to  a uniform  body  force,  between  parallel 
plates  with  stick  boundary  conditions.  If  the  plates  are  located  at  y = ±h  and  the  flow  is 
in  the  x direction,  then  according  to  linear  Navier-Stokes  hydrodynamics,  the  steady  state 
x-component  of  the  flow  velocity,  ux(y ),  satisfies 


ux(y) 


(2.1) 


where  Fe  is  the  driving  force  for  the  flow  (such  as  gravity),  n is  the  number  density,  and 
rj  is  the  shear  viscosity.  The  corresponding  temperature  profile  in  the  slit  is 


T(y)  = T0  - 


(■ nFe)2y4 
(12r,\) 


(2.2) 


where  To  is  the  midchannel  temperature  and  A is  the  thermal  conductivity  of  the  fluid. 
The  heat  generated  by  the  viscous  damping  is  dissipated  by  thermal  conduction  through 
the  walls  of  the  slit  [18]. 
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This  suggests  that  it  might  be  feasible  to  determine  both  the  thermal  conductivity  and  the 
shear  viscosity  of  a fluid  from  a single  nonequilibrium  molecular  dynamics  simulation  of 
Poiseuille  flow.  This  would  entail  determining  both  the  temperature  profile  and  the  flow 
velocity  profile.  Denis  Evans  and  co-workers  have  examined  the  issues  associated  with  the 
molecular  dynamics  simulation  of  Poiseuille  flow  in  a series  of  papers  [19,  20,  18,  21,  22,  23]. 
Their  objective  was  not  the  extraction  of  transport  coefficients  but  rather  understanding 
the  issues  involved  in  performing  steady-state  nonequilibrium  simulations.  Some  of  the 
salient  results  from  this  body  of  work  are  discussed  in  the  following  paragraphs. 

Their  simulations  used  a simple  arrangement  that  is  sketched  in  Fig.  2.1.  The  simulation 
cell  of  dimensions  LxxLyxLz  is  periodic  in  all  three  dimensions.  The  plates  are  represented 
by  three  layers  of  atoms  at  one  boundary  (dark  circles  in  the  figure)  that  are  tethered  to  fee 
sites  by  a harmonic  potential.  Note  that  only  one  set  of  layers  is  in  the  basic  simulation  cell. 
The  other  layer  is  just  the  periodic  image  of  the  first  layer.  The  wall  atoms  also  interact 
with  the  fluid  molecules  by  a specified  potential.  The  wall  atoms  are  thermostated  so 
that  the  heat  generated  in  the  fluid  is  removed  from  the  system.  The  fluid  atoms  are  not 
thermostated  and,  in  addition  to  intermolecular  interactions,  are  subject  to  a uniform,  time 
independent  driving  force  in  the  x direction.  The  use  of  a uniform  external  driving  force 
means  that  density  gradients  that  would  be  produced  by  a pressure  gradient  are  absent. 
Once  a steady  state  is  achieved,  the  velocity  profile,  ux(y)  and  the  temperature  profile  in 
the  fluid  are  determined  as  time  averages. 


O 

O 

o 

O O 

o 


• •••<> 
mu 


Fig.  2.1.  A sketch  of  the  simulation  scheme  used  by  Todd  and  Evans  [18].  The  dark  frame 
indicates  the  simulation  cell  (the  z-direction  is  normal  to  the  page)  and  the  dashed  lines  indicate 
the  effective  size  of  the  channel  (2h). 

The  initial  studies  used  a short  range  WCA-type  potential  to  describe  the  interactions 
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between  the  fluid  atoms  and  between  the  fluid  atoms  and  the  wall  atoms.  For  this  simple 
system,  it  was  necessary  to  postulate  the  existence  of  a strain  rate  term  in  the  heat  flux 
vector  in  order  to  accurately  describe  the  temperature  profiles  generated  in  the  simulations. 
In  addition  to  the  fourth  power  term  in  the  Navier-Stokes  solution,  this  coupling  generates 
a quadratic  term.  When  a molecular  fluid  (liquid  chlorine)  was  examined,  the  quadratic 
term  in  the  temperature  profile  was  not  observed  [23].  Also,  when  ethane  and  hexane  in 
fairly  narrow  channels  were  simulated,  the  quartic  form  described  the  temperature  profiles 
except  near  the  walls. 

A major  concern  in  many  of  these  studies  was  how  to  accurately  and  efficiently  generate 
the  temperature  profile.  The  kinetic  approach  uses  the  particle  velocity  minus  the  flow 
velocity  to  determine  the  local  temperature,  namely 

1 N 

kBT{y)  = £»«i(vi(i/)  - u(t/))2  (2.3) 

!# 

where  the  sum  is  over  the  molecules  at  a distance  y from  the  center  of  the  channel.  Some 
special  sampling  methods  were  developed  to  facilitate  evaluation  of  spatial  distributions 
without  recourse  to  the  “bins”  used  in  histogram  methods  [20].  The  kinetic  approach  is 
computationally  inefficient  for  nonequilibrium  simulations  as  it  requires  two  simulations, 
one  to  determine  the  flow  profile  u(y)  and  then  a second  simulation  to  determine  T(y)  or, 
if  the  distribution  is  generated  within  a single  simulation  by  taking  ratios  of  the  velocity 
at  a position  divided  by  the  density  at  that  position.  It  requires  long  runs  to  reduce  the 
statistical  scatter. 


An  alternative  method  for  determining  the  temperature  that  depends  only  on  the  configu- 
ration of  the  atoms  was  examined  in  some  detail.  The  theory  for  this  approach  is  described 
extensively  in  Ref.  [24]  and  is  based  on  the  thermodynamic  relation 


1 

k^T 


(-) 

\dEJv 


(2.4) 


The  details  are  not  pertinent  here,  so  we  only  consider  the  form  that  is  readily  implemented 
in  a simulation  without  the  need  for  multiple  simulations  or  extra  force  loops.  The  result 

is 


1 

k^f 


<-£v<-F*>/(£  f?>. 


(2.5) 


l 


l 


Another  concern  with  the  simulation  of  flow  in  a channel  is  how  the  width  of  the  channel 
is  reflected  in  departures  from  the  predictions  of  the  continuum  Navier-Stokes  equations 
[21].  Significant  departures  from  the  simple  hydrodynamics  predictions  are  found  when  the 
width  of  the  channel  (2h)  is  less  than  about  10  particle  diameters.  The  “stick”  boundary 
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condition  does  not  strictly  apply  when  the  wall  and  fluid  atoms  are  different  in  that  the 
fluid-fluid  and  fluid- wall  interactions  are  significantly  different.  The  lesson  appears  to  be 
that  behavior  near  the  walls  shows  significant  departures  from  the  predictions  of  Navier- 
Stokes  hydrodynamics. 

Finally,  the  examples  in  the  literature  are  for  short-range  potentials.  It  is  not  known 
what  the  result  would  be  if  longer  range  interactions  or  Coulomb  interactions  were  used. 
System  size  effects  would  be  of  particular  concern.  Coulomb  interactions  would  have 
additional  complications  since  periodicity  in  the  {/-direction  would  seem  to  be  unphysical. 
One  possible  approach,  that  would  not  be  too  expensive  computationally,  would  be  to  use 
the  Ewald  scheme  for  slab  geometry  that  has  a very  large  repeat  distance  in  the  {/-direction 
with  most  of  that  region  empty  [25].  Alternatively,  one  could  use  the  scheme  proposed  by 
Lekner  [26,  27]  although  it  is  computationally  quite  expensive. 

Implementation 

In  this  section,  we  describe  the  implemention  of  the  Poiseuille  Flow  simulation  and  also 
discuss  the  strengths  and  weaknesses  of  the  proposed  methods  based  on  our  experience. 
First  we  discuss  a method  that  follows  the  approach  of  Evans,  et  al.  Then  some  alternative 
methods  are  discussed.  The  examples  are  for  atomic  fluids  interacting  via  a Lennard-Jones 
12-6  potential.  Quantities  are  reported  in  units  of  the  Lennard-Jones  parameters  e and  a 
of  the  fluid  atoms  and  the  mass  of  the  particles  is  the  unit  of  mass  of  the  fluid  atoms. 

In  the  first  method  to  be  considered,  we  followed  the  approach  described  above  with 
one  modification.  Instead  of  placing  a wall  layer  at  one  end  of  the  simulation  cell  and 
considering  the  second  wall  .layer  to  be  the  periodic  image  of  the  first  layer,  we  considered 
having  two  explicit  walls  within  the  simulation  cell,  as  indicated  in  Fig.  2.2. 

Initially,  we  did  not  have  periodic  boundary  conditions  in  the  ^-direction.  However,  this 
was  not  fully  satisfactory  as  the  wall  layers  tended  to  be  pushed  out  by  the  fluid.  After 
some  experimentation,  we  settled  on  having  the  wall  layers  as  indicated  in  Fig.  2 but  with 
periodic  boundary  conditions  in  all  three  dimensions.  With  periodic  boundary  conditions, 
the  two  wall  layers  are  in  fact  one  larger  layer  with  two  surfaces  in  contact  with  the  fluid. 
This  eliminates  most  of  the  distortion  of  the  walls  by  the  fluid. 

The  code  written  to  simulate  Poiseuille  flow  contains  2000  Lennard-Jones  fluid  particles 
and  480  Lennard-Jones  wall  particles.  The  tethering  positions  for  the  wall  molecules  are 
specified  to  be  fee  sites  that  are  consistent  with  the  x-  and  z-dimensions  of  the  simulation 
cell.  The  L-J  a for  the  fluid- wall  interaction  is  set  to  1.2  to  strongly  inhibit  penetration 
of  the  wall  by  the  fluid  and  the  cr  for  the  wall- wall  interactions  is  set  to  2-1/6  so  that  the 
“equilibrium”  separation  matches  the  tethering  sites.  The  remainder  of  the  Lennard-Jones 
parameters  are  set  to  unity  for  all  interactions.  The  mass  of  all  particles  is  unity. 
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Fig.  2.2.  A sketch  of  the  simulation  scheme  used  here.  The  dark  frame  indicates  the  simulation 
cell  (the  z-direction  is  normal  to  the  page)  and  the  dashed  lines  indicate  the  effective  size  of  the 
channel  (2h).  Note  that  the  upper  wall  layers  are  within  the  simulation  cell,  unlike  Fig.  2.1. 

The  temperature  of  the  wall  particles  is  maintained  by  a Nose-Hoover  thermostat  [28].  The 
equations  of  motion  are  integrated  using  the  Beeman  algorithm  [2]  with  a time  step  of  O.Olr. 
Initially,  the  system  was  equilibrated  with  no  external  driving  force  for  a temperature  of 
1.0.  The  density,  n of  the  fluid  in  the  region  away  from  the  walls  settled  out  to  be  0.825. 

A series  of  runs  with  driving  forces  Fe  varying  from  0.01  to  0.05  were  implemented.  For  this 
system,  0.01  was  about  the  smallest  driving  force  that  resulted  in  a detectable  temperature 
profile.  The  kinetic  temperature  is  determined  at  each  time  step  using  the  local  flow 
velocity  (eq.  2.3)  at  that  instant.  Once  the  driving  force  was  applied,  experience  showed 
that  a run  of  2000r  (105  time  steps)  was  needed  to  obtain  a stationary  temperature  and 
flow  velocity  profile.  Then  a final  run  of  lOOOr  was  made  to  obtain  the  profiles  that  were 
fit  to  the  hydrodynamic  forms.  The  profile  information  is  collected  in  bins  of  width  lcr. 
An  example  of  the  profiles  is  displayed  in  Fig.  2.3. 

The  analysis  procedes  in  two  steps.  First,  the  flow  velocity  of  the  fluid  away  from  the  wall 
is  fit  to  a parabola  using  linear  least  squares, 

ux{yi)  = u0  4-  u2{yi  4-  yc)2  (2.6) 

where  yc.  the  ‘‘center”  of  the  fluid  layer  is  the  position  of  the  maximum  in  the  parabola 
and  yi  is  the  position  of  bin  i relative  to  yc.  The  centers  the  bins  in  the  simulation  cell  are 
Y = yi  + yc.  An  example  of  the  fit  is  shown  in  Fig.  2.4. 

This  simplifies  the  next  step  of  fitting  the  temperature  profile  T(yi)  to 

T(yi)  = r0  - T4yt,  (2.7) 
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F«  = 0.02,  Twa])  = 1.0 


Fig.  2.3.  The  profiles  for  the  density,  temperature,  flow  velocity  and  mean  square  variation  of  the 
flow  velocity  are  shown  for  the  case  with  Fx  = 0.02 


F*  = 0.02,  T^n  =1.0 


Fig.  2.4.  The  velocity  profile,  ux(Y),  is  shown  for  the  case  with  Fx  = 0.02.  The  filled  circles  are 
the  calculated  values  and  the  solid  curve  is  the  fitted  profile,  ux{Y)  — A + BY  4-  cY 2 . 

as  indicated  in  Fig.  2.5. 

The  transport  coefficients  are  then  determined  using 


—nFx 

2u2 

(2.8) 

( nFx )2 

(2.9) 

1277X4 
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and 


F,  = 0.02.  1.0 


Fig.  2.5.  The  temperature  profile,  T(y),  is  shown  for  the  case  with  Fx  — 0.02.  The  filled  circles 
are  the  calculated  values  and  the  solid  curve  is  the  fitted  profile,  T(y)  = Tq  — T^yA. 


The  results  for  several  different  values  of  the  driving  force  Fx  are  listed  in  Table  2. 

Table  2.  The  driving  force  Fx,  and  the  computed  values  for  the  shear  viscosity,  77  and  the  thermal 
conductivity,  A are  listed  for  the  case  with  the  density  n = 0.825  and  the  wall  temperature, 

Twall  = 1-0. 


Fx 

V 

A 

0.010 

2.57 

7.39 

0.015 

2.49 

7.34 

0.020 

2.50 

7.27 

0.025 

2.47 

7.44 

0.030 

2.42 

7.62 

0.040 

2.35 

8.00 

0.050 

2.32 

7.60 

These  estimates  for  77  and  A are  consistent  with  published  values  of  these  transport  coeffi- 
cients obtained  using  the  time  correlation  method  of  statistical  mechanics  [4,  7]. 

In  order  to  get  some  idea  for  the  uncertainties  associated  with  these  estimates  for  77  and 
A,  a set  of  ten  runs  of  lOOOr  duration  were  made  for  a driving  force,  Fx  = 0.025,  for  the 
same  state  as  identified  in  Table  2.  The  viscosity  values  ranged  from  2.43  to  2.49  with 
an  average  and  uncertainty  (the  variance)  of  2.45T0.018.  The  themal  conductivity  values 
ranged  from  7.11  to  8.32  with  an  an  average  and  variance  of  7.46T0.36. 

In  order  to  obtain  some  idea  of  how  long  a simulation  should  run,  the  state  identified  in 
Table  2 was  run  for  200r.  The  cumulative  time  averages  for  77,  A,  and  n are  shown  in 
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Fig.  2-6.  The  convergence  for  77  is  quite  good  while  the  convergence  for  A is  not  convincing 
for  lOOOr. 


t/l 


Fig.  2.6.  The  cumulative  values  for  77,  A,  and  the  number  density  n are  shown  a simulation  of 
2000t.  The  wall  temperature  is  Twau  = 1,  and  Fx  = 0.025  The  degree  of  convergence  for  77  is 
very  good,  while  that  for  A is  not  as  definite. 


Some  further  tests  were  made  by  generating  time  averages  of  lOr  over  an  interval  of  5000r. 
The  shear  viscosity  values  were  closely  grouped  around  the  average  while  the  thermal 
conductivity  values  were  broadly  distributed  over  the  interval  from  5 to  12.  From  this  we 
conclude  that  this  method  can  generate  the  shear  viscosity  from  a lOOOr  run  while  much 
longer  runs  are  needed  to  produce  thermal  conductivity  estimates  with  reasonably  small 
uncertainties.  This  is  consistent  with  the  observation  that  long  runs  are  needed  to  reduce 
statistical  scatter  when  the  kinetic  definition  of  local  temperature  is  used. 

In  order  to  reduce  the  scatter  in  the  thermal  conductivity  values,  another  set  of  ten  sim- 
ulations was  performed  with  the  local  kinetic  temperature  being  determined  using  the 
accumulated  flow  profile,  u(y),  rather  than  the  instanteneous  value.  The  duration  of  each 
of  these  simulations  was  2000r  as  this  was  found  to  be  necessary  for  approximate  con- 
vergence of  the  thermal  conductivity.  The  averages  and  variances  for  the  ten  runs  are 
77  = 2.38±0.01,  A = 7.34±0.40,  and  n - 0.832±0.00. 
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Table  3.  The  results  for  the  shear  viscosity,  the  thermal  conductivity  and  the  number  density  are 
listed  for  the  ten  simulations  where  the  cumulative  flow  velocity  was  used  to  determine  the  kinetic 
temperature. 


V 

A 

n 

2.35 

7.09 

0.832 

2.36 

7.39 

0.832 

2.38 

7.50 

0.832 

2.39 

7.70 

0.832 

2.39 

7.25 

0.832 

2.40 

8.19 

0.832 

2.42 

6.76 

0.832 

2.35 

7.33 

0.832 

2.37 

7.39 

0.832 

2.38 

6.76 

0.832 

Alternative  boundary  regions. 

Two  variants  of  the  boundary  region  were  tested.  The  approach  of  Evans  is  rather  inflexible 
in  that  the  wall  atoms  are  tethered  to  specific  sites.  There  is  some  ambiguity  about  where 
the  fluid  meets  the  wall.  Also,  the  process  of  changing  the  density  of  the  fluid  by  changing 
the  dimensions  of  the  simulation  cell  involves  some  additional  work.  For  these  reasons, 
some  alternatives  were  considered. 

Suppose  the  regions  at  the  ends  of  the  cell  are  considered  to  be  the  wall  region,  as  shown  in 
Fig.  2.2.  Molecules  in  the  wall  region  are  fluid  molecules  and  are  free  to  move  in  and  out  of 
the  wall  region.  If  a molecule  reaches  either  end  of  the  cell,  it  is  diffusely  reflected  back  into 
the  fluid  with  momenta  selected  from  a Boltzmann  distribution  with  the  wall  temperature 
imposed.  This  proved  to  be  rather  inefficient  in  establishing  a good  temperature  profile 
and  was  discontinued. 

Instead,  the  approach  that  has  been  selected  treats  the  two  bins  at  the  ends  of  the  cell  as 
the  “wall” . Periodic  boundary  conditions  are  imposed  in  all  directions.  At  each  time  step, 
the  molecules  in  the  wall  region  are  thermostated  to  the  wall  temperature  by  rescaling  the 
momenta,  and  a zero  net  flow  condition  in  the  direction  of  the  applied  force  that  is  imposed 
on  the  molecules.  In  other  respects,  this  approach  is  the  same  as  the  other  methods. 
This  leads  to  coarse  grained  temperature  and  flow  profiles  that  consistent  are  with  those 
predicted  by  linear  hydrodynamics.  The  shear  viscosity  and  thermal  conductivity  are 
estimated  by  fitting  quadratic  and  quartic  profiles. 

An  extensive  set  of  runs  to  implement  this  method  has  been  made  for  a Lennard-Jones 
fluid  with  a mean  density  of  0.85  and  average  temperatures  ranging  from  0.9  to  2.75.  The 
results  are  displayed  in  Fig.  2.7  for  systems  that  were  stabilized  for  2000r  followed  by  a 
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production  run  of  an  additional  2000r.  The  viscosities  are  in  reasonable  agreement  with 
those  listed  in  Table  1.  For  temperatures  greater  than  2,  the  thermal  conductivities  are 
larger  than  those  in  Table  1.  Also,  the  scatter  in  the  thermal  conductivity  values  is  large. 


Poiseuille  Flow  for  transport  coefficients 
L-J  fluid,  n=0.85 


Fig.  2.7.  The  shear  viscosity  and  thermal  conductivity  for  a range  of  states  for  the  Lennard- Jones 
fluid  are  shown  here.  The  circles  are  for  Fx  =0.3,  the  other  symbols  are  for  F^  =0.4.  Varying  the 
number  of  layers  in  the  system  from  19  to  35  has  minimal  effect  on  the  results. 

3.  Momentum  Inpulse  Relaxation 

This  approach  to  the  determination  of  the  shear  viscosity  uses  the  decay  in  time  of  a care- 
fully designed  transient  momentum  gradient  and  extracts  the  shear  viscosity  by  matching 
the  simulated  momentum  gradient  to  the  hydrodynamic  solution  for  the  same  problem 
[29].  The  essence  of  the  approach  is  the  following. 

Suppose  the  flow  velocity  is  ux(y , t)  given  that  the  initial  value  is  ux(y , 0)  = ao  exp (—b0y2). 
If  there  is  no  pressure  gradient  or  external  force  present,  the  flow  velocity  satisfies 

dux{y,t)  d2ux(y,t) 

dt  V dy2  ( ' 

The  solution  is 

Ux(y,t)=  + 1/^}  1/2exp[-fr0j/2/(l  + t/t0)]  (3.2) 
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where  to  = l/(4i/6o),  and  v is  the  kinematic  viscosity,  the  shear  viscosity  divided  by  the 
density. 

While  this  appears  to  be  a straightforward  simulation,  there  is  a problem  with  the  periodic 
boundary  conditions  in  the  ^/-direction.  If  the  usual  condition  of  not  changing  the  y- 
component  of  the  momentum  when  a molecule  leaves  the  simulation  cell  through  a plane 
perpendicular  to  the  y- axis  and  its  periodic  image  enters  the  cell  through  the  opposite  side 
of  the  cell,  the  flow  velocity  will  not  decay  to  zero  at  long  times.  To  take  care  of  this,  a 
special  procedure  is  introduced  to  permit  the  flow  velocity  to  decay  as  if  the  system  were 
infinite.  This  is  done  by  requiring  the  x-component  of  the  velocity  of  the  particle  entering 
the  simulation  cell  to  be  consistent  with  the  current  estimate  of  the  Gaussian  velocity 
profile. 

The  process  of  determining  the  Gaussian  profile  at  a time  t is  achieved  by  partitionaing 
the  simulation  cell  into  bins  parallel  to  the  x-z  plane.  The  flow  velocity  of  each  bin  is  just 
the  average  x-component  of  the  velocity  of  the  particles  in  each  bin  at  that  time.  The 
resulting  velocity  profile  is  then  fitted  to  a two  parameter  Gaussian  function, 

ux(t)  = up(t)  exp(-y2  / a2  (t)) . (3.3) 


The  fit  parameters,  up(t)  and  cr2(£),  are  then  used  to  determine  by  extrapolation  the  x- 
component  of  the  velocity  of  the  image  molecule,  uf'(i),  just  before  it  enters  the  simulation 
cell  as 

uf  (t)  = up(t)  exp (~y2{t  - A t)/a2{t))  (3.4) 

where  y(t  — At)  is  the  ^/-coordinate  of  the  image  molecule  at  the  time  step  before  it  enters 
the  cell.  The  decay  of  the  imposed  flow  velocity  is  realized  by  assigning  u^(t)  to  the 
x-component  of  the  velocity  of  the  particle  entering  the  cell. 


This  fitting  procedure  is  also  used  to  determine  the  shear  viscosity.  Both  up(t)  and  cr2{t) 
can  be  used  to  estimate  the  characteristic  time  to  and  therefore  the  viscosity.  The  peak 
velocity,  up(t)  has  the  form 


A 

{l  + Bt)1/2 


(3.5) 


so  that 


pB 
4 bn 


(3.6) 


Also,  the  variance  cr2(t)  is  a linear  function  of  time  and  may  be  used  to  provide  another 
estimate  of  to- 


In  the  reported  tests,  the  cell  had  dimensions  L/2  x L x L/2  and  19  bins  were  used.  The 
dimensions  were  such  that  on  the  order  of  100  molecules  were  in  each  bin.  The  simulations 
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were  performed  by  equilibrating  the  fluid  without  any  imposed  flow  field,  and  then  running 
the  system  with  the  initially  imposed  flow  for  about  10  ps.  This  sequence  was  repeated 
20  times  with  separate  initial  conditions.  Results  are  available  for  the  Lennard-Jones 
fluid  and  for  “butane”.  The  results  for  the  shear  viscosity  compare  favorably  with  those 
obtained  both  with  Green-Kubo  time  correlations  and  with  nonequilibrium  applied  shear 
simulations.  The  required  run  times  for  this  method  are  at  least  and  order  of  magnitude 
smaller  that  for  the  other  methods. 

Some  possible  concerns  are  about  the  uniformity  of  the  density  (not  discussed)  and  about 
appropriate  parameters  for  the  imposed  initial  conditions  (extensively  discussed  below). 
There  has  not  been  any  proposal  for  a similar  approach  to  determining  the  thermal  con- 
ductivity of  a fluid,  although  the  decay  of  a temperature  pulse  has  been  used  in  a molecular 
dynamics  simulation  to  estimate  the  thermal  conductivity  of  a superlattice  crystal  [30].  It 
would  require  a careful  analysis  of  thermal  expansion  effects  to  implement  such  an  approach 
for  a liquid. 

Implementation 

The  Momentum  Impulse  Relaxation  scheme  has  been  implemented  as  follows.  First,  a set 
of  coordinates  for  2000  molecules  at  the  desired  temperature  and  density  are  generated 
with  a uniform  momentum  distribution.  Then  a momentum  pulse  with  amplitude  a( 0) 
and  width  6(0)  are  imposed  and  the  system  evolves  for  lOr.  The  evolution  of  a(t ) and  b(t) 
are  determined  using  19  bins  uniformly  distributed  in  the  simulation  cell.  This  sequence 
is  repeated  20  times  and  the  average  of  a(t)  is  used  to  determine  r\  using  eq.  3.5.  Some 
results  for  n = 0.85  and  three  values  of  a(0)  (0.3,  0.4,  and  0.5)  and  6(0)  = 0.25  are  shown 
in  Fig.  3.1.  If  a(0)  were  smaller  than  0.3,  it  was  difficult  to  extract  a useful  signal.  These 
values  for  the  viscosity  are  consistent  with  those  listed  in  Table  1. 

While  it  is  straightforward  to  set  up  the  momentum  pulse  in  an  initially  “equilibrium” 
fluid  state,  I have  encountered  serious  difficulties  in  realizing  results  in  the  form  of  eq.  3.2 
over  extended  times.  The  task  is  to  determine  the  amplitude  and  width  of  the  momentum 
pulse  by  fitting  the  calculated  profile  to  a Gaussian  form.  This  is  a nonlinear  fitting  task, 
but  since  the  changes  from  one  time  step  to  the  next  are  small,  this  is  not  a difficult  thing 
to  do  using  a steepest  descent  method.  For  relatively  short  times,  the  fitted  amplitude  and 
width  conform  the  the  scaled  forms.  However,  for  times  longer  than  one  time  unit  or  so, 
there  are  serious,  systematic  departures  from  the  scaled  form. 

My  conjecture  is  that  the  suggested  way  of  dealing  the  the  “leakage”  of  momentum  out  of 
the  simulation  cell  is  not  adequate.  Actually  only  a few  particles  cross  the  cell  boundary 
at  any  one  time  so  the  flux  out  of  the  simulation  cell  is  small.  Since  this  method,  while 
requiring  relatively  short  overall  simulation  times  is  potentially  useful,  it  is  not  at  this  time 
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Fig.  3.1.  The  estimated  values  of  the  shear  viscosity  and  indicated  for  three  values  of  the  initial 
amplitude  of  the  momentum  pulse.  The  circles  are  for  a(0)  = 0.3,  squares  are  for  a(0)  = 0.4,  and 
the  diamonds  are  for  a(0)  = 0.5.  The  lines  connecting  the  points  are  to  guide  the  eye. 


sufficiently  robust  to  warrant  further  investigation  here. 

4.  Reversed  perturbation  nonequilibrium  MD 

Another  approach  to  obtaining  transport  coefficients  from  molecular  simulations  is  based 
on  the  following  observations.  A transport  coefficient  is  the  proportionality  coefficient 
between  a current  and  a force.  For  example, 


(4.1) 


where  A is  the  thermal  conductivity,  dT/dz  is  the  temperature  gradient  that  drives  the 
current,  and  J is  the  heat  current.  In  the  “usual”  molecular  dynamics  approaches  to  deter- 
mining a transport  coefficient,  the  driving  force  is  specified  and  the  current  is  calculated. 
In  the  reversed  perturbation  approach,  the  current  is  established  by  some  scheme  and  the 
resulting  gradient  is  determined  [31,  32,  33,  34].  In  these  papers,  the  methodology  for 
determining  thermal  conductivity  and  shear  viscosity  are  examined.  Some  initial  applica- 
tions for  the  determination  of  the  thermal  conductivity  for  a variety  of  systems  have  been 
reported  [35,  36,  37,  38,  39,  40,  41].  Extensions  of  the  method  to  diffusion  in  mixtures  has 
also  been  described  [42,  43]. 

The  essential  feature  of  this  approach  is  to  generate  a known  current  that  is  stationary  in 
time.  This  is  done  by  separating  the  simulation  cell  of  dimensions  Lx  x Ly  x Lz  into  layers 
as  indicated  in  Fig.  4.1.  The  central  layer  and  the  layer  at  one  end  of  the  cell  are  special. 
At  some  interval  rtrans , a property  of  a molecule  in  each  of  these  layers  is  interchanged 
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such  that  energy  and  momentum  are  conserved  and  that  a gradient  between  the  layers  is 
established.  Because  of  periodic  boundary  conditions  in  all  three  dimensions,  a current 
flows  from  the  central  layer  to  the  layer  at  one  end  and  also  to  the  periodic  image  of  the 
end  layer.  Once  a stationary  current  is  established,  one  then  determines  the  corresponding 
gradient  by  sampling  the  temperature  or  flow  velocity  in  each  of  the  layers. 


Fig.  4.1.  A sketch  of  the  simulation  scheme.  The  dark  frame  indicates  the  simulation  cell  (the 
x-direction  is  normal  to  the  page).  The  dashed  lines  indicate  the  boundaries  of  the  layers.  The 
cross  hatched  layers  are  the  end  layer  and  its  periodic  image  and  the  tilted  line  layer  is  the  central 
layer.  Note  that  the  total  number  of  layers  is  even  so  the  the  distance  between  the  “walls”  and 
the  central  cell  is  the  same  in  both  directions. 


First  let  us  consider  the  way  a heat  current  is  established.  Let  the  central  layer  be  the  hot 
one  and  the  other  layer  be  the  cold  one.  At  time  intervals  of  Ttrans,  the  velocity  vectors 
of  the  molecule  with  the  largest  kinetic  energy  of  molecules  in  the  cold  layer  and  the 
molecule  with  the  lowest  kinetic  energy  in  the  hot  layer  are  exchanged.  This  conserves  the 
total  momentum  and  energy  of  the  fluid  and  eventually  results  in  a temperature  gradient 
between  the  central  layer  and  the  two  crosshatched  layers  shown  in  Fig.  3.1.  The  heat 
current  is  then  the  difference  in  the  kinetic  energy  transfered  between  the  hot  and  cold 
layers  averaged  over  time  and  over  the  area  of  the  layer.  If  v\  and  v\  are  the  squares  of 
the  velocities  of  the  hot  and  cold  molecules,  then  the  time  averaged  heat  current,  J,  is 


J = 


__  Etrans/ers”^-^)/2 


2 //  Zy  j*  Jj  -j 


(4.2) 


where  t is  the  averaging  time  and  transfers  is  the  set  of  energy  transfers.  If  an  exchange 
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of  kinetic  energy  occurs  once  every  rtrans  then  t/rtrans  is  the  number  of  transfer  events. 

The  temperature  gradient  is  obtained  from  the  temperature  profile  of  the  layers.  The 
temperature  of  layer  k,  Tfc,  is  2/3  times  the  average  kinetic  energy  of  the  Nk  molecules  in 
layer  k divided  by  Boltzmann’s  constant  ks- 


Tk 


3 Nkks 


(4.3) 


Determining  the  average  temperature  of  the  layers  is  less  time  consuming  that  evaluat- 
ing the  time  correlation  function  of  the  properly  defined  heat  current  that  goes  into  the 
Green-Kubo  expression  for  the  thermal  conductivity.  The  time  consuming  part  is  the 
establishment  of  a stationary  temperature  gradient. 

A similar  scheme  is  used  to  establish  a momentum  flux.  Instead  of  picking  the  hottest  and 
coldest  molecules  in  the  central  and  end  layers,  one  picks  the  molecule  with  the  largest 
^/-component  of  velocity  to  the  right  in  the  central  layer  and  the  molecule  with  the  largest 
^component  of  velocity  in  the  opposite  direction  in  the  end  layer  and  interchanges  the 
velocity  vectors.  Again  the  total  momentum  and  energy  are  conserved.  The  result  is  a 
gradient  of  the  z-direction  of  the  ^/-component  of  the  layer  velocity.  The  layer  velocity,  uk 
is  just  the  average  velocity  of  the  molecules  in  the  k- th  layer.  The  momentum  transfer 
for  a transfer  event,  A py,  is  just  the  difference  in  the  ^components  of  momentum  of  the 
selected  molecules  and  the  shear  viscosity,  77,  is  then 


_ Ytransfer  APy  (4.4} 

2 tLxLydu/dz 

and  the  shear  gradient,  du/dz , is  the  slope  of  the  uy(k)  vs  z(k)  curve. 

Some  limits  on  the  applicability  of  this  approach  have  been  noted  [44,  40].  If  the  switching 
time  rtrans  is  too  long,  then  the  gradients  do  not  develop  as  intended  and  the  simple 
connection  between  the  flux  and  the  temperature  or  velocity  profile  is  lost. 

Implementation 

This  approach  requires  a separate  simulation  for  the  determination  of  the  shear  viscosity 
and  the  thermal  conductivity.  First,  we  examine  the  results  for  the  thermal  conductivity, 
A,  for  the  state  with  T = 1.0  and  the  density  n = 0.887.  The  results  for  a set  of  ten 
(10)  sequential  simulations  of  duration  lOOOr  (105  time  steps)  are  listed  in  Table  4.  The 
temperature  profile  for  one  of  these  runs  is  displayed  in  Fig.  4.2  to  illustrate  the  sort  of 
profiles  involved  in  determining  dT(y)/dy. 

The  corresponding  results  for  the  shear  viscosity  are  found  in  Fig.  4.3  and  Table  5. 
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A separate  pair  of  runs,  one  each  for  77  and  for  A,  of  106At  were  also  made,  the  results 
are:  A P = 0.0598,  77  = 3.75  and  3.76;  Jq  = 0.0737,  A = 8.76  and  8.79.  The  results  of  the 
longer  runs  are  consistent  with  those  listed  in  Tables  4 and  5. 

As  a further  check  on  this  method,  a series  of  runs  were  made  for  the  states  listed  in 
Table  1.1.  The  results  for  n = 0.85  are  displayed  in  Fig.  4.4. 


y/o 

Fig.  4.2.  The  temperature  profile,  T(y ),  is  shown  for  a case  with  switching  interval  =50  At.  The 
solid  circles  are  the  computed  profile  and  the  solid  lines  are  the  least-squares  fit  to  the  profile. 


Table  4.  The  results  for  the  thermal  conductivity  are  listed  for  10  runs,  each  of  lOOOr  duration; 
Jq  = 0.0735T0.004  and  A = 8.53T0.26.  The  two  values  of  the  thermal  conductivity  listed  for 
each  flux  are  determined  from  the  two  parts  of  the  profile  shown  in  Fig.  4.2.  The  switching  interval 
is  50  At. 


Jq 

A 

A 

0.0737 

8.40 

8.44 

0.0732 

8.05 

8.16 

0.0738 

8.99 

8.97 

0.0741 

8.56 

8.54 

0.0735 

8.54 

8.60 

0.0730 

8.34 

8.20 

0.0740 

8.78 

8.69 

0.0730 

8.35 

8.38 

0.0731 

8.47 

8.44 

0.0738 

8.80 

8.97 
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y/a 

Fig.  4.3.  The  momentum  profile,  u(y),  is  shown  for  a case  with  switching  interval  = 50 At.  The 
solid  circles  are  the  computed  profile  and  the  solid  lines  are  the  least-squares  fit  to  the  profile. 


Table  5.  The  results  for  the  shear  viscosity  are  listed  for  10  separate  runs, each  of  lOOOr  duration; 
AP  = 0.0599±0.0002  and  rjs  = 3.80T0.13.  The  two  values  of  the  shear  viscosity  listed  for  each 
flux  are  determined  from  the  two  parts  of  the  profile  shown  in  Fig.  4.3.  The  switching  interval  is 
1~trans  — 5(3  At. 


A P 

V 

V 

0.0601 

4.05 

4.01 

0.0598 

3.88 

3.95 

0.0596 

3.55 

3.61 

0.0599 

3.76 

3.66 

0.0598 

3.85 

3.79 

0.0601 

3.78 

3.87 

0.0600 

3.93 

3.90 

0.0599 

3.70 

3.69 

0.0601 

3.85 

3.74 

0.0598 

3.72 

3.78 

The  interval  between  swap  events,  rtrans  is  an  important  parameter  that  determines  the 
magnitude  of  the  temperature  or  momentum  gradient.  Some  tests  of  the  sensitivity  of 
the  results  have  been  made  for  the  T = 3.564,  n = 0.65  state  of  Table  1.  The  results, 
that  were  generated  by  Jim  Rainwater  [45],  are  listed  in  Table  6.  These  results  indicate 
that  0.25r  < rtransfer  < 1.0 r is  an  acceptable  range  of  switching  times.  The  Green-Kubo 
estimates  for  this  state  are  A = 5.4±0.3  and  rjs  = 1.15T0.03. 
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9.5 
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Fig.  4.4.  The  results  of  this  method  (black  squares)  are  compared  with  the  Green-Kubo  results 


(open  circles)  of  Table  1.1  for  the  n = 0.85  states. 


Table  6.  The  results  for  the  shear  viscosity  and  thermal  conductivity  as  a function  of  the  switching 
time  Ttrans  are  listed  for  the  T = 3.564,  n = 0.65  state  of  Table  1. 


Ttrans  / T 

A 

r] 

0.06 

5.06 

0.12 

5.12 

1.13 

0.25 

5.20 

1.16 

0.40 

5.19 

0.50 

5.17 

1.16 

1.00 

5.21 

1.19 

2.00 

5.11 

1.22 

4.00 

5.11 

1.21 

5.  Discussion 

First,  we  consider  the  relative  wall  clock/CPU  time  needed  to  obtain  “satisfactory”  es- 
timates for  the  shear  viscosity  and  the  thermal  conductivity  using  the  various  methods 
examined.  For  purposes  of  reference,  we  note  that  a Green-Kubo  simulation  of  duration 
104t  generates  estimates  with  approximately  5%  uncertainty  in  the  shear  viscosity  and 
about  10%  uncertainty  in  the  thermal  conductivity. 

As  indicated  by  the  values  in  Table  3 and  Figure  2.6,  the  results  for  Poiseuille  flow  are 
that  a simulation  of  103r  will  produce  shear  viscosity  estimates  with  an  uncertainty  of  less 
than  5%.  Even  allowing  for  the  larger  system  size  (a  factor  of  5),  this  requires  less  CPU 
time  than  the  Green-Kubo  method.  That  is  not  the  case  for  the  thermal  conductivity  as 
a simulation  nearly  20  times  longer  is  needed  to  achieve  a 10%  uncertainty  level. 
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For  the  Reversed  Pertubation  Nonequilibrium  Molecular  Dynamics  method,  the  results  in 
Tables  4 and  5 suggest  that  the  5%  uncertainty  for  the  shear  viscosity  and  10%  uncertainty 
for  the  thermal  conductivity  can  be  obtained  with  103r  simulations.  Now  allowing  for  the 
need  to  establish  the  gradients  and  the  larger  system  size  involved,  the  individual  transport 
coefficients  can  be  obtained  in  a time  comparable  to  that  required  for  the  Green-Kubo 
method.  If  both  coefficients  are  needed,  the  Green-Kubo  method  would  be  a bit  faster. 

The  Momentum  Impulse  Relaxation  method  is  considerably  faster  than  any  of  the  other 
methods.  However,  until  the  problem  of  departures  from  the  scaling  of  the  solutions  from 
the  hydrodynamic  predictions  are  resolved,  this  method  remains  a research  topic. 

There  are  some  “loose  ends”  with  the  nonequilibrium  simulation  methods.  These  are 
the  “optimal”  switching  rate  for  the  Reversed  Perturbation  method  and  the  appropriate 
magnitude  of  the  driving  force  for  Poiseuille  flow.  Another,  potentially  more  significant 
factor  is  the  minimum  system  size  needed  to  obtain  proper  profiles.  The  simulations 
discussed  here  contained  over  2000  atoms.  If  this  could  be  reduced  to  1000,  the  advantages 
of  the  nonequilibrium  methods  over  the  Green-Kubo  method  would  be  significant. 

Another  factor  to  be  considered  concerns  problems  associated  with  these  methods  when 
applied  to  molecular  fluids  as  opposed  to  the  monatomic  Lennard-Jones  fluid.  If  we  restrict 
our  attention  to  rigid  molecules,  then  there  do  not  appear  to  be  serious  problems  provided 
the  momenta  and  forces  are  center-of-mass  values  and  molecular  currents  are  used  [46,  47]. 
Further  discussion  of  this  topic  lies  outside  the  scope  of  this  report. 

Finally,  we  note  that  cases  where  the  viscosity  of  the  fluid  is  sufficiently  large  that  vis- 
coelastic effects  are  present  have  not  been  examined.  For  those  situations,  it  is  known 
that  the  correlation  times  for  the  shear  viscosity  correlation  functions  appearing  in  eq.  1.4 
become  much  longer  than  occurs  for  any  of  the  cases  examined  here.  In  that  situation,  the 
Green-Kubo  method  will  require  much  longer  simulations  in  order  to  obtain  statistically 
independent  samples.  It  is  not  known  how  the  other  methods  will  behave  for  such  systems. 
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